##################################################################################
##################################################################################
####### Values and Political Predispositions in the Age of Polarization: #########
#### Unpacking the Relationship between Partisanship and Ideology, 1988-2012 #####
##################################################################################
##################################################################################
  
####
## The following code is used to 
## both examine changes in the reliability
## of the value orientations scale over time,
## and conduct an item analysis for purposes
## of checking that the item response
## functions for each value item
## are, in fact, roughly montonic
## (as assumed by the summated rating model).
####


.libPaths("/Library/Frameworks/R.framework/Versions/3.0/Resources/library")

library(stats)
library(car)
library(lattice)
library(MiscPsycho)
library(MASS)
library(foreign)
library(mokken)

setwd("/Users/adamenders/Dropbox/Values/")


#################################################################################################

### Examine changes in value scale reliability (Cronbach's Alpha) over time ###

val.items <- read.dta("")
head(val.items)
ncol(val.items)

## Partition value items by year

sub88 <- subset(val.items[,2:11], val.items$Time==0)
sub92 <- subset(val.items[,2:11], val.items$Time==1)
sub96 <- subset(val.items[,2:11], val.items$Time==2)
sub00 <- subset(val.items[,2:11], val.items$Time==3)
sub04 <- subset(val.items[,2:11], val.items$Time==4)
sub08 <- subset(val.items[,2:11], val.items$Time==5)
sub12 <- subset(val.items[,2:11], val.items$Time==6)


## Compute the reliability of the scale for each year

check.reliability(val.items[,2:11])$alpha 


alpha <- c(check.reliability(sub88)$alpha,
           check.reliability(sub92)$alpha,
           check.reliability(sub96)$alpha,
           check.reliability(sub00)$alpha,
           check.reliability(sub04)$alpha,
           check.reliability(sub08)$alpha,
           check.reliability(sub12)$alpha
)

alpha.data <- as.data.frame(alpha)

year <- c(1988, 1992, 1996, 2000, 2004, 2008, 2012)

rownames(alpha.data) <- year


## Dotplot of Cronbach's Alpha for Value Scale Over Time

dotplot(year ~ alpha, data = alpha.data,
        aspect = 1.5,
        xlab = "Cronbach's Alpha",
        ylab = "",
        xlim = c(.62, .82),
        col = 1,
        lyt = 2,
        scales=list(
          y=list(labels = rownames(alpha.data))),
        panel=function(x, y){
          panel.dotplot(x, y, col = "black", lty = 2)
        }
)


#################################################################################################

### Conduct item analysis by examing item response functions for value items ###

val.items <- read.dta("")
head(val.items)

valuesdata <- val.items[,2:11]

sumrate <- as.data.frame(scale(valuesdata))

sd(sumrate)

# Use "alpha" function to create an object containing parts of a reliability analysis.
reliab.sumrate <- alpha(sumrate)
reliab.sumrate$coefficients

sumrate <- as.data.frame(scale(val.items[,2:11]))

sd(sumrate)

reliab.sumrate <- alpha(sumrate)
reliab.sumrate$coefficients

####
## Use "apply" function to create summated rating scale of 
## issue responses, defined as the mean of each person's 
## responses across the ten items
####

scale <- apply(val.items[,2:11], MARGIN = 1, FUN = mean)

## Summary statistics on the scale

summary(scale)
sd(scale)

## Histogram of original scale

histogram(~scale, aspect = .75,
          xlab = "Values Summated Rating Scale",
          col = "grey")

####
## Perform an item analysis on the scale
## Print alpha coefficients for restscores, removing each variable, 
## in turn, from overall ten-item scale.
####

reliab.sumrate$condAlpha

# Calculate restscores and plot against corresponding items

restscores <- matrix(nrow = nrow(val.items[,2:11]), ncol = ncol(val.items[,2:11]))

names <- rep("rest", times = ncol(valuesdata))

for (j in 1:ncol(valuesdata)) {
  restscores[, j] <- apply(valuesdata[, -j], MARGIN = 1, FUN = mean)
  names[j] <- paste(names[j], colnames(valuesdata)[j], sep = ".")
}

colnames(restscores) <- names

set.seed(1357)
xyplot(valuesdata[,1] ~ restscores[,1],
       aspect = 1,
       main = "A. Equalopp",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,2] ~ restscores[,2],
       aspect = 1,
       main = "B. Equalrights",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,3] ~ restscores[,3],
       aspect = 1,
       main = "C. Equalchance",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,4] ~ restscores[,4],
       aspect = 1,
       main = "D. Lessequal",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,5] ~ restscores[,5],
       aspect = 1,
       main = "E. Unequal",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,6] ~ restscores[,6],
       aspect = 1,
       main = "F. Fewer",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,7] ~ restscores[,7],
       aspect = 1,
       main = "G. Changing",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,8] ~ restscores[,8],
       aspect = 1,
       main = "H. Lifestyles",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,9] ~ restscores[,9],
       aspect = 1,
       main = "I. Standards",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)


set.seed(1357)
xyplot(valuesdata[,10] ~ restscores[,10],
       aspect = 1,
       main = "J. Family",
       ylab = "Response Categories",
       xlab = "Rest Scores",
       panel = function (x, y) {
       panel.xyplot(jitter(x), jitter(y), col = "gray")
       panel.spline(x, y, col = "black")
       }
)

